﻿using System;

namespace NCM_MSTest.Alg
{
    /// <summary>
    /// 最小二乘拟合
    /// </summary>
    public class LeastSquaresFit
    {
        /// <summary>
        /// 直线最小二乘拟合 y=a+bx
        /// 输入：实验数据x，y
        /// 输出：拟合参数最佳估计值a、b，标准偏差sigma,相关系数r
        /// </summary>
        /// <param name="x">实验数据x</param>
        /// <param name="y">实验数据y</param>
        /// <param name="yest">直线最小二乘拟合得到的函数y的值</param>
        /// <param name="a">拟合参数最佳估计值</param>
        /// <param name="b">拟合参数最佳估计值</param>
        /// <param name="sigma">标准偏差</param>
        /// <param name="r">相关系数</param>
        public static void LinearFit(double[] x, double[] y
            , ref double[] yest
            , ref double a, ref double b
            , ref double sigma, ref double r)
        {
            double sx = 0, sy = 0, xy = 0, xx = 0, yy = 0, q = 0;
            int n = x.Length;
            for (int i = 0; i < x.Length; i++)
            {
                sx = sx + x[i];
                sy = sy + y[i];
                xy = xy + x[i] * y[i];
                xx = xx + x[i] * x[i];
                yy = yy + y[i] * y[i];
            }
            a = (sy * xx - sx * xy) / (n * xx - sx * sx);
            b = (n * xy - sx * sy) / (n * xx - sx * sx);
            r = (xy - sx * sy / n) / Math.Sqrt((xx - sx * sx / n) * (yy - sy * sy / n));
            for (int i = 0; i < x.Length; i++)
            {
                yest[i] = a + b * x[i];
                q = q + Math.Pow(y[i] - yest[i], 2);
            }
            sigma = Math.Sqrt(q / (n - 1));
        }

        #region 可转换为直线关系的模型

        /// <summary>
        /// 幂函数：y=d*x^b
        /// </summary>
        /// <param name="x"></param>
        /// <param name="y"></param>
        /// <param name="yest"></param>
        /// <param name="a"></param>
        /// <param name="b"></param>
        /// <param name="sigma"></param>
        /// <param name="r"></param>
        public static void LinearFit_PowerFunction(double[] x, double[] y
            , ref double[] yest
            , ref double a, ref double b
            , ref double sigma, ref double r)
        {
            double[] _x = new double[x.Length];
            double[] _y = new double[x.Length];
            for (int i = 0; i < x.Length; i++)
            {
                _x[i] = Math.Log(x[i]);
                _y[i] = Math.Log(y[i]);
            }
            LinearFit(_x, _y, ref yest, ref a, ref b, ref sigma, ref r);
            double d = Math.Pow(10, a);
            a = d;
            for (int i = 0; i < x.Length; i++)
            {
                yest[i] = d * Math.Pow(x[i], b);
            }
        }

        /// <summary>
        /// 指数函数：y=d*e^{b*x}
        /// </summary>
        /// <param name="x"></param>
        /// <param name="y"></param>
        /// <param name="yest"></param>
        /// <param name="a"></param>
        /// <param name="b"></param>
        /// <param name="sigma"></param>
        /// <param name="r"></param>
        public static void LinearFit_ExponentialFunction(double[] x, double[] y
            , ref double[] yest
            , ref double a, ref double b
            , ref double sigma, ref double r)
        {
            double[] _x = new double[x.Length];
            double[] _y = new double[x.Length];
            for (int i = 0; i < x.Length; i++)
            {
                _x[i] = x[i];// Math.Log(x[i], Math.E);
                _y[i] = Math.Log(y[i], Math.E);
            }
            LinearFit(_x, _y, ref yest, ref a, ref b, ref sigma, ref r);
            double d = Math.Pow(Math.E, a);
            a = d;
            for (int i = 0; i < x.Length; i++)
            {
                yest[i] = d * Math.Pow(Math.E, b * x[i]);
            }
        }

        /// <summary>
        /// 对数函数：y=a+b*logx
        /// </summary>
        /// <param name="x"></param>
        /// <param name="y"></param>
        /// <param name="yest"></param>
        /// <param name="a"></param>
        /// <param name="b"></param>
        /// <param name="sigma"></param>
        /// <param name="r"></param>
        public static void LinearFit_LogarithmicFunction(double[] x, double[] y
            , ref double[] yest
            , ref double a, ref double b
            , ref double sigma, ref double r)
        {
            double[] _x = new double[x.Length];
            double[] _y = new double[x.Length];
            for (int i = 0; i < x.Length; i++)
            {
                _x[i] = Math.Log(x[i]);
                _y[i] = y[i];// Math.Log(y[i], Math.E);
            }
            LinearFit(_x, _y, ref yest, ref a, ref b, ref sigma, ref r);
            for (int i = 0; i < x.Length; i++)
            {
                yest[i] = a + b * Math.Log(x[i]);
            }
        }

        /// <summary>
        /// 双曲函数：1/y = a+ b* 1/x
        /// </summary>
        /// <param name="x"></param>
        /// <param name="y"></param>
        /// <param name="yest"></param>
        /// <param name="a"></param>
        /// <param name="b"></param>
        /// <param name="sigma"></param>
        /// <param name="r"></param>
        public static void LinearFit_HyperbolicFunction(double[] x, double[] y
            , ref double[] yest
            , ref double a, ref double b
            , ref double sigma, ref double r)
        {
            double[] _x = new double[x.Length];
            double[] _y = new double[x.Length];
            for (int i = 0; i < x.Length; i++)
            {
                _x[i] = 1 / x[i];
                _y[i] = 1 / y[i];
            }
            LinearFit(_x, _y, ref yest, ref a, ref b, ref sigma, ref r);
            double d = Math.Pow(Math.E, a);
            a = d;
            for (int i = 0; i < x.Length; i++)
            {
                yest[i] = 1 / (a + b * 1 / x[i]);
            }
        }

        #endregion 可转换为直线关系的模型

        public static double Determinant(double[,] A)
        {
            if (A.GetLength(0) != A.GetLength(1))
                throw new ArgumentException();
            return 0;
        }

        public static void CramersRule(double[,] A, double[] B, ref double[] k)
        {

        }

        /// <summary>
        /// n阶多项式最小二乘拟合
        /// </summary>
        /// <param name="n">多项式阶数</param>
        /// <param name="x">实验数据x</param>
        /// <param name="y">实验数据y</param>
        /// <param name="yest">拟合得到的函数y的值</param>
        /// <param name="a">拟合参数最佳估计值</param>
        /// <param name="sigma">标准偏差</param>
        /// <param name="r">相关系数</param>
        public static void PolyFit(int n, double[] x, double[] y
            , ref double[] yest
            , ref double[] a
            , ref double sigma, ref double r)
        {
            a = new double[n];
            int N = x.Length;
            double[,] A = new double[n, n];
            double[] B = new double[n];
            //A x a = B 根据克莱默法则计算系数数组a
            A[0, 0] = N;
            for (int i = 0; i < n; i++)
            {
                for (int j = 0; j < n; j++)
                {
                    if (i == 0 && j == 0)
                        continue;
                    A[i, j] = 0;
                    for (int k = 0; k < N; k++)
                    {
                        //计算矩阵A
                        A[i, j] = A[i, j] + Math.Pow(x[k], j);
                    }
                }
                //计算矩阵B
                B[i] = B[i] + y[i] * Math.Pow(x[i], i);
            }
        }
    }
}